'''
Author: chyang0822 270917365@qq.com
Date: 2025-03-31 09:06:58
LastEditors: chyang0822 270917365@qq.com
LastEditTime: 2025-03-31 11:02:06
FilePath: /Project-DMAI/guide/pybind11/myadd.py
Description: 

Copyright (c) 2025 by ${git_name_email}, All Rights Reserved. 
'''
# myadd.py
import numpy as np
import scipy
import skimage
import warnings
import vtk
def read_image(image_path):
    """读取医学图像并返回VTK图像数据及参数"""
    reader = vtk.vtkNIFTIImageReader()
    reader.SetFileName(image_path)
    reader.Update()
    vtk_image_data = reader.GetOutput()
    # 获取QForm矩阵
    qform_matrix = reader.GetQFormMatrix()
    if qform_matrix is None:
        raise ValueError("无法获取QForm矩阵。请检查NIfTI文件的有效性。")
    # 提取原点、方向矩阵和spacing
    origin = np.array([qform_matrix.GetElement(i, 3) for i in range(3)])
    direction_matrix = np.array([[qform_matrix.GetElement(i, j) for j in range(3)] for i in range(3)])
    spacing = np.array(vtk_image_data.GetSpacing())
    
    return vtk_image_data, origin, direction_matrix, spacing

def convert_to_numpy_and_filter(vtk_image_data, sigma):
    """将VTK图像转换为numpy数组并应用高斯滤波"""
    from vtkmodules.util import numpy_support
    from scipy.ndimage import gaussian_filter
    
    dims = vtk_image_data.GetDimensions()
    vtk_array = vtk_image_data.GetPointData().GetScalars()
    data = numpy_support.vtk_to_numpy(vtk_array)
    data = data.reshape(dims, order='A')
    
    # 高斯平滑处理
    return gaussian_filter(data, sigma=sigma)

def calculate_threshold(data, p1):
    """计算图像分割的灰度阈值"""
    hist, bin_edges = np.histogram(data.flatten(), bins='auto')
    cumulative_hist = np.cumsum(hist)
    total_pixels = cumulative_hist[-1]
    threshold_index = np.where(cumulative_hist / total_pixels > p1)[0][0]
    return bin_edges[threshold_index]

def segment_and_label(data, threshold):
    """图像分割和连通域标记"""
    from scipy.ndimage import label as nd_label
    from skimage.measure import regionprops
    
    segmented_data = np.where(data > threshold, 1, 0)
    labeled_array, num_features = nd_label(segmented_data)
    regions = regionprops(labeled_array)
    
    return labeled_array, regions

def extract_spheres(regions, origin, direction_matrix, spacing, 
                     min_radius_mm, max_radius_mm, min_sphericity):
    """从区域中提取符合条件的球体"""
    detected_spheres = []
    
    for region in regions:
        # 获取质心坐标
        ijk_centroid = np.array(region.centroid)
        world_centroid = direction_matrix @ (ijk_centroid * spacing) + origin

        # 计算半径
        coords = region.coords
        distances = np.linalg.norm((coords - ijk_centroid) * spacing, axis=1)
        fitted_radius_mm = np.max(distances)

        # 计算圆度
        convex_hull_area = region.convex_area * (spacing[0] * spacing[1])
        actual_area = region.area * (spacing[0] * spacing[1])
        circularity = actual_area / convex_hull_area if convex_hull_area > 0 else 0

        # 筛选条件检查
        if min_radius_mm <= fitted_radius_mm <= max_radius_mm and circularity >= min_sphericity:
            print(f"连通域 {region.label}:")
            print(f"  - 半径: {fitted_radius_mm:.2f} mm")
            print(f"  - 圆度: {circularity:.2f}")
            
            # 创建球体对象并添加到列表
            # sphere = DetectedSphere(world_centroid, fitted_radius_mm, circularity)
            # detected_spheres.append(sphere)
    
    return detected_spheres

def sphere_detection(image_path, p1=0.99997, sigma=0.5, min_radius_mm=1, max_radius_mm=3.0, min_sphericity=0.7):
    import warnings
    warnings.filterwarnings('ignore')
    
    # 读取图像
    vtk_image_data, origin, direction_matrix, spacing = read_image(image_path)
    
    # 将VTK图像转换为numpy数组并进行高斯滤波
    data_array = convert_to_numpy_and_filter(vtk_image_data, sigma)
    
    # 计算阈值并进行分割
    threshold = calculate_threshold(data_array, p1)
    print(f'计算的灰度阈值 T: {threshold:.2f}')
    
    # 图像分割和连通域分析
    labeled_array, regions = segment_and_label(data_array, threshold)
    
    # 提取符合条件的球体
    detected_spheres = extract_spheres(regions, origin, direction_matrix, spacing, 
                                                min_radius_mm, max_radius_mm, min_sphericity)
    
def SupportVersion(a, b):
    print("Now is in python module--------")
    print("numpumpy version : ",np.__version__)
    print("scipy version : ",scipy.__version__)
    print("skimage version : ",skimage.__version__)
    print("vtk version : ",vtk.__version__)
    print("{} + {} = {}".format(a, b, a+b))
    print("Now is in python module--------------;-")
    return a + b